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Abstract 

The T-matrix method is widely used for the calculation of scattering by particles of sizes on the 
order of the illuminating wavelength. Although the extended boundary condition method (EBCM) is 
the most commonly used technique for calculating the T-matrix, a variety of methods can be used. 

We consider some general principles of calculating T-matrices, and apply the point-matching 
method to calculate the T-matrix for particles devoid of symmetry. This method avoids the time- 
consuming surface integrals required by the EBCM. 
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1 The T-matrix method 

The T-matrix method in wave scattering involves writing the relationship between the wave incident 
upon a scatterer, expanded in terms of orthogonal eigenfunctions, 

00 fincl 

Umc = Y.""^n , (1) 

n 

where a„ are the expansion coefficients for the incident wave, and the scattered wave, also expanded in 
terms of orthogonal eigenfunctions, 

oo 

Uscat = tPm ' (2) 

k 

where p^ are the expansion coefficients for the scattered wave, is written as a simple matrix equation 

oo 

Pk = Y, '^knan (3) 
n 

or, in more concise notation, 

P = TA (4) 

where Tj-„ are the elements of the T-matrix. The T-matrix method can be used for scalar waves or vector 
waves in a variety of geometries, with the only restrictions being that the geometry of the problem 
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permits expansion of the waves as discrete series in terms of orthogonal eigenfunctions, that the response 
of the scatterer to the incident wave is linear, and that the expansion series for the waves can be truncated 
at a finite number of terms. In general, one calculates the T-matrix, although it is conceivable that it might 
be measured experimentally. 

The T-matrix depends only on the particle — its composition, size, shape, and orientation — and is 
independent of the incident field. This means that for any particular particle, the T-matrix only needs to 
be calculated once, and can then be used for repeated calculations. This is a significant advantage over 
many other methods of calculating scattering where the entire calculation needs to be repeated 1 1 1. Some 
cases provide even more efficiency: if the waves are expanded in spherical functions, the averaging of 
scattering over various orientations of the particle compared to the direction of the incident wave can be 
performed analytically |2|. 

In the spherical geometry of elastic light scattering by a particle contained entirely within some ra- 
dius ro, the eigenfunction expansions of the fields are made in terms of vector spherical wavefunctions 

(vswFs) (HElElliEIEIlzl: 

Mi'/)(lcr) = N„hi''^\kr)Cnm{9,cl)) (5) 

N.(^e?(^0-^^^^^]B.^ (6) 

where h„ ' (kr) are spherical Hankel functions of the first and second kind, N„ = 1/ \/n{n + 1) are 
normalisation constants, and Bnm{d,4>) = r VY;f(g,0), Cn m{d,(p) = V x (rY/f (0, c^))), and Pnm{d,4)) = 
fy^(0, (p) are the vector spherical harmonics I1I2I3I4I5I6I7I . and 0) are normalised scalar spherical 
harmonics. The usual polar spherical coordinates are used, where 9 is the co-latitude measured from the 
+z axis, and cp is the azimuth, measured from the +x axis towards the +y axis. 

M^m and N^^J, are outward-propagating TE and TM multipole fields, while M^^J, and n|,^2 are the 
corresponding inward-propagating multipole fields. Since these wavefunctions are purely incoming 
and purely outgoing, each has a singularity at the origin. Since fields that are free of singularities are of 
interest, it is useful to define the singularity-free regular vector spherical wavefunctions: 

RgM„„,(fcr) = l[Ml'i(/cr)+Ml'i(/cr)], (7) 
RgN„„,(fcr) = i[Nll,l(fcr) + Nlil(fcr)]. (8) 

Since the spherical Bessel functions {kr) = ^ (/zi^^ (kr) + hl^^ {kr)), the regular VSWFs are identical to the 
incoming and outgoing VSWFs except for the replacement of the spherical Hankel functions by spherical 
Bessel functions. 

Since the incident field, in the absence of a scatterer, is singularity-free, the expansion 

oo f! 

Einc(r)=I I a\?J,RgM„^{kr)+b\?J,RgN„^{kr) (9) 

n=l m=—n 



is generally used for the incident field. Alternatively, the purely incoming part of the incident field can 
be used: 

, ,,(2)^yr(2)/Z,.^^^,{2)xT(2)/ 

n=l m= 

In both cases, the scattered field is 



Einc(r)=£ £ fll'iMl'i(fcr)+bii|Nl'i(fcr). (10) 



Escat(r) = £ £ pl!]Ml^i(/cr) +4mNl^i(fcr). (11) 

(!=1 m=— J! 
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The two sets of expansion coefficients for the incident /incoming field are related, since fl„m = 2flJ,m 

(3) (2) 

and = 2b However, the scattered /outgoing field expansion coefficients will differ, as will the 
T-matrix. Using the regular expansion, the T-matrix in the absence of a scatterer is a zero matrix, while 
using the incoming field expansion, the no-scatterer T-matrix is the identity matrix. The two expan- 
sions are essentially the same — the only difference is that the incident wave in the incident/ scattered 
wave expansion includes part of the outgoing wave. T-matrices for the two expansions only differ by 
the identity matrix, so t(^"/°"*) = 2T(™'^/*'^^*) + I. The incident /scattered formulation is much more 
commonly used; the incoming/ outgoing formulation gives simpler results for the transport of momen- 
tum and angular momentum (that is, optical force and torque) by the field. It should be note that for 
plane wave illumination, for which the VSWF expansion is non-terminating, the incident/ scattered for- 
mulation gives a scattered wave expansion that converges over all space, while the incoming/ outgoing 
expansion, strictly used, would give an non-terminating, non-convergent outgoing field expansion. For 
focussed beam illumination with a finite VSWF expansion, the incoming /outgoing expansion directly 
gives the total outgoing field that would be experimentally measured. Since conversion from one formu- 
lation to the other is simple, either can be readily used for calculation of fields, forces, scattering matrices, 
or for orientation averaging. 

In practice, the field expansions and the T-matrix are terminated at some n = Nmax- For the case of 
a scatterer that is contained within a radius vq, Nmax ~ kro is usually adequate, but Nmax = kro + 3\/kro 
is advisable if higher accuracy is needed Q. Although we assume in this paper (as is usually the case) 
that the incident and scattered wave expansions are terminated at the same Nmax (giving a square T- 
matrix), this is not necessary. It should be noted that convergence of the expansion of the incident 
field is not a necessary condition for the T-matrix method to be useful — indeed, for the most common 
application, scattering of an incident plane wave, the incident field expansion does not converge over 
all space. However, it does converge within the radius — which is the part of the field that can affect 
the scattering particle — and therefore, the field expansions and the T-matrix can be truncated at a finite 

Nmax- 

For the case of plane wave scattering, the plane wave expansion formula is useful: 

Unn, = 47ri"N„C:,„ • Eo, b„m = 47ri"-iN„B*„, ■ Eq. (12) 

The main case of interest for non-plane wave incident illumination is that of focussed beams. A 
variety of methods can be used, such as plane wave expansion L8J, the localised approximation L9 ^10llTTl 
IT2l . or the point-matching method fl^l . 

The only remaining requirement is that the T-matrix be calculated. This requires essentially a com- 
plete calculation of the scattering properties of the scatterer. This is almost universally done using the ex- 
tended boundary condition method (EBCM), originally developed by Waterman 1 3 1, which is so strongly 
linked with the T-matrix method that the terms "EBCM" and "T-matrix method" are often used inter- 
changeably. In the next section, we consider some general principles involved in the calculation of the T- 
matrix, and show that an alternative method — column-by-column calculation using the point-matching 
method (PMM) — is computationally feasible and simply implemented for homogeneous isotropic parti- 
cles devoid of symmetry. 

Lastly, before we continue to consider calculation of T-matrices in more detail, we can note that while 
the incident and scattered fields are usually expanded in terms of VSWFs, other sets of eigenfunctions, 
such as cylindrical wavefunctions (for scatterers of infinite length in one dimension), or a Floquet ex- 
pansion (planar periodic scatterers), are more appropriate for other geometries. There is no requirement 
that the modes into which the incident and scattered fields are expanded be the same, or even similar. In 
all of these cases, the T-matrix method remains applicable. 
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2 Calculating the T-matrix 



If the field expansions and T-matrix are truncated at some Nmax, there are Nj = 2Nmax(Nmax + 2) 
expansion coefficients for each of the incident and scattered fields, and the T-matrix is Nt x Nj- Since 
Nmax is proportional to the radius enclosing the particle, ro, the number of expansion coefficients is 
proportional to Tq, and the number of elements in the T-matrix is proportional to r^. This can be used to 
obtain an estimate of the scaling of computational time for different methods of calculation. 

2.1 The extended boundary condition method 

In principle, any method of calculating scattering by the particle can be used to calculate the T-matrix. 
However, the method of choice is almost universally the EBCM |ll|2||3j|4|. In the EBCM, the internal 
field within the particle is expanded in terms of regular VSWFs. Therefore, the method is restricted to 
homogeneous and isotropic particles. Rather than considering the coupling of the incident and scat- 
tered fields directly, the coupling between the incident and internal (the RgQ matrix), and scattered and 
internal fields (the Q matrix) is calculated, and the T-matrix found from these (T = — RgQQ^^). The 
RgQ and Q matrices are the same size as the T-matrix, with 0{N^^^) elements. The elements of these 
matrices are found by integrating over the surface of the scatterer, an operation requiring 0{N^^^) time 
per element, so the calculation of the RgQ and Q matrices is expected to require 0{N^^^) computational 
time. The actual calculation of the T-matrix, if direct inversion is naively used, takes 0{N^^^) time. In 
practice, the calculation of the RgQ and Q matrices dominates the computational time 1141 . 

From this, it can be seen that the EBCM can be expected to be very slow for large particles. However, 
most applications of the EBCM have been for the special case of scattering particles rotationally symmet- 
ric about the z axis. In this case, scattered modes m|J]j, and n|j!]j, only couple to incident modes RgM^^ 
and RgN„„^ if m' = m, greatly reducing the number of matrix elements that need to be calculated, and 
the surface integral over the particle surface reduces to a one-dimensional integral over 9, since the az- 
imuthal integration over (p can be simply done analytically |4|. This results in a great improvement in 
performance, and, in terms of computational time, EBCM is clearly the method of choice for axisym- 
metric particles. Numerical problems do occur when the scatterer is highly non-spherical. The discrete 
sources method is designed to overcome these problems fl5|. For the even more symmetric case of a 
spherical scatterer, the scattered and incident modes only couple if n' = n and m' = m, the T-matrix 
becomes diagonal, and all of the integrals can be performed analytically, and Mie's solution to scattering 
by a sphere 1 16 1 is simply obtained. 

In a similar manner, scatterers with point-group rotational symmetry allow significant improvement 
of the computational time required through exploitation of the symmetry Il^'l/'^ISV 

Methods have also been developed to calculate T-matrices for clusters of particles and for layered 
particles C). 

The efficiency of the EBCM for the calculation of the T-matrix is such that alternative methods 
need only be considered if the EBCM is inapplicable (such as when the particle in inhomogeneous or 
anisotropic), numerical difficulties are encountered using the EBCM (such as for extremely non-spherical 
particles), or if the scattering particle has no symmetries that can be used to optimise the computation of 
the T-matrix. 

2.2 Methods other than the EBCM 

Methods other than the EBCM can be used to calculate the T-matrix. In general, one would calculate the 
scattered field, given a particular incident field. The most direct way in which to use this to produce a 
T-matrix is to solve the scattering problem when the incident field is equal to a single spherical mode — 
that is, a single VSWF such as Einc(r) = RgM^-^{kr), Einc(r) = RgN;^;^(fcr), Einc(r) = RgM2;^(A:r), etc, 
and repeat this for all VSWFs that need to be considered (up to n = Nmax)- The expansion coefficients 
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for the scattered field can be found in each case, if necessary, by using the orthogonal eigenfunction 
transform (the generalised Fourier transform), and each scattering calculation gives a single column of 
the T-matrix. 

Therefore, the calculation of a T-matrix requires that 2Nmax(Nmax + 2) separate scattering problems 
are solved. The provides a criterion for deciding whether it is desirable to calculate a T-matrix: if more 
than 2Nmax(Nmax + 2) scattering calculations will be performed, then it is more efficient to calculate 
the T-matrix and use this for the repeated calculations than it is to use the original scattering method 
repeatedly. Repeated calculations are expected if orientation averaging is to be carried out, or if inho- 
mogeneous illumination is to be considered, such as, for example, scattering by focussed beams, where 
there are generally 6 degrees of freedom, namely the three-dimensional position of the scatterer within 
the beam, and the three-dimensional orientation of the scatterer. Even if only a modest number of points 
are considered along each degree of freedom, the total number of scattering calculations required rapidly 
becomes very large, and even if the T-matrix takes many hours to calculate, the total time saved by doing 
so can make an otherwise computationally infeasible problem tractable. 

Volume methods are of interest, since they can readily be used for inhomogeneous or anisotropic 
particles. The two most likely candidates are the finite-difference time-domain method (FDTD) 1 1 , 19J 
and the discrete dipole approximation (DDA). In FDTD, the Maxwell equations are discretised in space 
and time, and, beginning from a known initial state, the electric and magnetic fields at each spatial grid 
point are calculated for successive steps in time. The number of grid points required is O(N^ax) fo^" 
three-dimensional scattering, and O(Nmax) time steps required, so FDTD solutions scale as 0{N^^^). 
Therefore, calculation of the T-matrix using FDTD should scale as 0{N^^^), which is the same scaling 
as the EBCM. However, the grid required must be closely spaced compared to the wavelength, and 
the space outside the scatterer must also be discretised, making FDTD substantially slower than EBCM, 
especially for smaller particles. However, FDTD is an extremely general technique, and has potential as 
a method for the calculation of T-matrices. 

We should add that there is an additional consideration that makes FDTD potentially attractive as a 
method for calculating the T-matrix: FDTD does not assume that the incident wave is monochromatic. 
Consider the case when the illumination is a brief pulse with a Gaussian envelope. The frequency spec- 
trum of the incident wave is Gaussian, and the scattering of a range of frequencies can be found by taking 
the Fourier transform of the scattered field |20|. Even if we are not interested in other than monochro- 
matic illumination, we will frequently be interested in scattering by size distributions of particles. Since 
varying the frequency for a particular particle is equivalent to varying the size of the particle for a fixed 
incident frequency, the T-matrices for a range of particle sizes can be calculated simultaneously. 

The other major volume method for computational scattering, the discrete dipole approximation 
(DDA), also known as the coupled-dipole method, has been recently applied to the calculation of the T- 
matrix by Mackowski |21|, who obtained good results, with reasonable computational efficiency using 
a moment method to solve the DDA system of equations. DDA lacks the main disadvantages of FDTD, 
namely the need to discretise space outside the particle, and the need to implement suitable boundary 
conditions to prevent non-physical reflections from the boundary of the computational domain. Mack- 
owski's method scales as 0{N'^^^) for large Nmax- There is no need to discuss his method in detail here, 
and the interested reader is referred to his recent description of the method |21 1. 

Finally, we consider the point-matching method. Like the T-matrix method and the EBCM, the point- 
matching method involves expansion of fields in terms of VSWFs. In the point-matching method, the 
internal field within the scatterer and the scattered field are expanded as series of VSWFs, and the inci- 
dent, internal, and scattered fields are matched at points on the particle surface, using the usual bound- 
ary condition of continuity of tangential components of the electric and magnetic fields. This gives a 
system of equations from which the unknown expansion coefficients of the internal and scattered fields 
can be found. Typically, enough points are used for matching the fields so as to give an overdetermined 
system of equations, which is then solved in a least-squares sense. Solving the system of 0{N^^^) un- 
knowns using direct matrix inversion can be expected to be an O(N^ax) problem, with the result that the 
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total computational time is ©(N^ax)- Iri practice, faster methods can be used, and our results indicate a 
performance of about O(N^ax) fo^" implementation. 

The point-matching method is an attractive candidate since a T-matrix implementation will generally 
include routines to calculate VSWFs, making the implementation of a point-matching T-matrix calculator 
simple. The only further requirement is a routine for solving overdetermined linear systems. Since the 
scattered field is calculated in terms of VSWFs, the conversion of the results of a single PMM calculation 
to a T-matrix column is trivial. 

Naturally, multiple expansion methods (the generalised multipole technique, or the multiple multi- 
pole method) can be used. Since multipole methods exist for anisotropic media 1 22 1, the method can be 
used for anisotropic scatterers. 

Our implementation of the point-matching method, and its performance, is discussed in the next 
section. 



3 Point-matching method 

Our implementation of the PMM T-matrix calculation uses an incoming/outgoing field expansion (equa- 
tions ilOll and lITTl ), rather than the usual incident/ scattered wave expansion (equations © and iTTll '). and 
the internal field is expanded in terms of regular VSWFs: 

Mr^(/cr)+fo„,„Nr„U/cr), (13) 

(1=1 m= — n 

Nmax n 

Escat(r) = £ £ p„„,MU(/:r) +(j„„,NU(/:r), (14) 

n=l m— — n 

Nmax n 

Eint(r) = £ I c„,„RgM„,„(/cr)+4mRgN„„(/cr). (15) 

)i=l ni— — n 

We use this particular expansion since we are interested in calculating optical forces and torques within 
optical traps |23 .24| and this results in simpler expressions for these quantities. 

We considered a single scatterer, centred on the origin, contained entirely within a radius ro, and with 
a surface specified by a function of angle: 

r = r(0,(/)) (16) 
The boundary conditions — matching the tangential fields on the surface of the scatterer — are 

n X (Einc(r) + Escat(r)) = n x Eint(r), (17) 
nx (Hinc(r) + Hscat(r)) = n x Hint(r), (18) 

where fi is a unit vector normal to the surface of the particle. 

The magnetic fields are given by expansions similar to those for the electric fields: 

1 ^"""^ " (2) (2) 

Hinc(r) = T X L anm^nm{kr) +bnmM^{kr), (19) 

"^medium fj—i m=—n 

Hscat(r) = T £ £ p„mN^m(fcr) -F(?„mMU(fcr), (20) 

"^medium n—1 m=—n 
-I Nmax n 

Hint(r) = £ £ c„,„RgN„„,(fcr)+4mRgM„,(/cr). (21) 

'^particle „=i m=-n 

where ^medium arid k^anicle are the wavenumbers of the field in the surrounding medium and inside the 
particle, respectively. 
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^''Omax 


Time 


1 


0.033 


0.041 s 


2 


0.21 


0.16 s 


3 


0.55 


0.85 s 


4 


1.00 


7.00 s 


5 


1.54 


30.3 s 


6 


2.14 


1.86 mm 


7 


2.78 


4.95 min 


8 


3.46 


12.2 min 


9 


4.17 


26.8 min 


10 


4.90 


56.3 min 


11 


5.66 


1.91 h 


12 


6.42 


3.53 h 


13 


7.21 


6.35 h 



Table 1: Computation times for calculating T-matrices. The calculations were carried out in MATLAB 
on a 1.5 GHz PC. The maximum size parameter krQ for which the truncation is expected to always be 
well-convergent is shown. Reasonable convergence can also be expected for size parameters kr^ ^ Nmax 

There are 4]Vmax(N'inax + 2) unknown variables — the expansion coefficients c„m, dnm, Pnm, and q„m- 
Since the fields are vector fields, each point gives multiple equations — four independent equations per 
point. We generate a grid of 2Nmax(Nmax + 2) points with equal angular spacings in each of the 9 and 
(p directions, giving 8Nmax(Ninax + 2) independent equations. Equal angle spaced points are used for 
simplicity, although points uniformly distributed about a sphere would be better 1251 . 

The values of the VSWFs at these points on the particle surface are calculated, and used in the 
column-by-column calculation of the T-matrix. 

The computation time (which is independent of the particle shape, depending only on the containing 
radius tq) is shown in tabled The calculations were carried out in MATLAB 1261 on a 1.5 GHz PC. 
The times taken are reasonable in comparison to computation times for EBCM for particles with no 
symmetry |14|. 

Results of a sample calculation are shown in figure 121 where the diagonal scattering matrix elements 
are shown, calculated using the PMM T-matrbc. The scattering matrix elements Sn and S22 are shown 
for scattering in two different planes; the effect of non-axisymmetry is evident. 

The accuracy and validity of the PMM-calculated T-matrix will be essentially the same as the ac- 
curacy and validity of the point-matching algorithm used in the calculation. Thus, a detailed analysis 
of our simple proof-of-principle implementation serves little purpose. It is obviously useful to use the 
best, sufficiently fast, point-matching code available. In view of the mathematical similarity between the 
T-matrix method and the point-matching method, it should be a simple task to adapt any PMM code to 
the task of generating T-matrix columns. 

4 Conclusions 

The point-matching method is suitable for the calculation of the T-matrix for particles with no symmetry, 
provided that the particles are not too large. The method has the advantage of being extremely simple to 
implement within a general T-matrix package, since most of the required routines will be shared with the 
existing T-matrix code. This results from the mathematical formalisms of the T-matrix method and the 
point-matching method being essentially the same. Any point-matching algorithm can be used, with 
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max 



Figure 1: Computation times for calculating T-matrices. The calculations were carried out in MATLAB 
on a 1.5 GHz PC. The time taken scales as 0{N^l^) for Nmax > 2. 

multiple expansion origins, automatic convergence checks, and so on. Since the PMM uses the same 
field expansions as the EBCM, the same numerical difficulties are to be expected for scatterers with large 
aspect ratios; in such cases, multiple expansion origin algorithms will be necessary. The accuracy of the 
PMM T-matrix will be the same as the PMM which is used to calculate it. Naturally, the usual conditions 
of applicability of the PMM, such as the validity of the Rayleigh hypothesis, need to be considered. 

The PMM explicitly depends on the Rayleigh hypothesis — the assumption that the fields can be rep- 
resented by the expansions llT3l — lITSt over all space rather than just outside and inside spherical surfaces 
circumscribing and inscribing the surface of the scatterer. The validity of this assumption for arbitrary 
scatterers is unknown. However, the use of an overdetermined system of equation may well extend the 
method somewhat beyond the strict range of applicability of the Rayleigh hypothesis by providing a 
least squares approximation of the fields between the circumscribing and inscribing surfaces where the 
VSWF expansions might be non-convergent. One advantage of relying on the Rayleigh hypothesis is 
that the fields are given everywhere, including the fields internal to the scatterer (a T^^"'' -matrix can be 
used to relate the internal and incident fields). This applies generally to methods that make use of the 
Rayleigh hypothesis, such as the generalised separation of variables method 1271 . In contrast to this, the 
EBCM, which avoids the Rayleigh hypothesis, gives the tangential surface fields on the surface of the 
scatterer, rather than the internal fields. 

The point-matching method lacks the generality of DDA and FDTD. In this respect, the recent dis- 
crete dipole moment method T-matrix calculations by Mackowski 1211 are particularly promising. 

Lastly, we note again that FDTD may prove to be a useful method for T-matrix calculation since it 
can be used to calculate T-matrices simultaneously for a range of particle sizes. 
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